      SUBROUTINE E02ADF(M,KPLUS1,NROWS,X,Y,W,WORK1,WORK2,A,S,IFAIL)
C
C     NAG LIBRARY SUBROUTINE  E02ADF
C
C     E02ADF  COMPUTES WEIGHTED LEAST-SQUARES POLYNOMIAL
C     APPROXIMATIONS TO AN ARBITRARY SET OF DATA POINTS.
C
C     FORSYTHE-CLENSHAW METHOD WITH MODIFICATIONS DUE TO
C     REINSCH AND GENTLEMAN.
C
C     USES NAG LIBRARY ROUTINE  P01AAF.
C     USES BASIC EXTERNAL FUNCTION  SQRT.
C
C     STARTED - 1973.
C     COMPLETED - 1976.
C     AUTHOR - MGC AND JGH.
C
C     WORK1  AND  WORK2  ARE WORKSPACE AREAS.
C     WORK1(1, R)  CONTAINS THE VALUE OF THE  R TH  WEIGHTED
C     RESIDUAL FOR THE CURRENT DEGREE  I.
C     WORK1(2, R)  CONTAINS THE VALUE OF  X(R)  TRANSFORMED
C     TO THE RANGE  -1  TO  +1.
C     WORK1(3, R)  CONTAINS THE WEIGHTED VALUE OF THE CURRENT
C     ORTHOGONAL POLYNOMIAL (OF DEGREE  I)  AT THE  R TH
C     DATA POINT.
C     WORK2(1, J)  CONTAINS THE COEFFICIENT OF THE CHEBYSHEV
C     POLYNOMIAL OF DEGREE  J - 1  IN THE CHEBYSHEV-SERIES
C     REPRESENTATION OF THE CURRENT ORTHOGONAL POLYNOMIAL
C     (OF DEGREE  I).
C     WORK2(2, J)  CONTAINS THE COEFFICIENT OF THE CHEBYSHEV
C     POLYNOMIAL OF DEGREE  J - 1  IN THE CHEBYSHEV-SERIES
C     REPRESENTATION OF THE PREVIOUS ORTHOGONAL POLYNOMIAL
C     (OF DEGREE  I - 1).
C
C     NAG COPYRIGHT 1975
C     MARK 5 RELEASE
C     MARK 6 REVISED  IER-84
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C
C     CHECK THAT THE VALUES OF  M  AND  KPLUS1  ARE REASONABLE
C
C     .. Parameters ..
      CHARACTER*6       SRNAME
      PARAMETER         (SRNAME='E02ADF')
C     .. Scalar Arguments ..
      INTEGER           IFAIL, KPLUS1, M, NROWS
C     .. Array Arguments ..
      DOUBLE PRECISION  A(NROWS,KPLUS1), S(KPLUS1), W(M), WORK1(3,M),
     *                  WORK2(2,KPLUS1), X(M), Y(M)
C     .. Local Scalars ..
      DOUBLE PRECISION  ALPIP1, BETAI, BJ, BJP1, BJP2, CI, D, DF, DI,
     *                  DIM1, DJ, EPSR, FACTOR, PIJ, SIGMAI, WRPR,
     *                  WRPRSQ, X1, XCAPR, XM
      INTEGER           I, IERROR, IPLUS1, IPLUS2, J, JPLUS1, JPLUS2,
     *                  JREV, K, MDIST, R
C     .. Local Arrays ..
      CHARACTER*1       P01REC(1)
C     .. External Functions ..
      INTEGER           P01ABF
      EXTERNAL          P01ABF
C     .. Intrinsic Functions ..
      INTRINSIC         SQRT
C     .. Executable Statements ..
      IERROR = 4
      IF (KPLUS1.LT.1 .OR. M.LT.KPLUS1) GO TO 380
      K = KPLUS1 - 1
C
C     TEST THE VALIDITY OF THE DATA.
C
C     CHECK THAT THE WEIGHTS ARE STRICTLY POSITIVE.
C
      IERROR = 1
      DO 20 R = 1, M
         IF (W(R).LE.0.0D0) GO TO 380
   20 CONTINUE
C
C     CHECK THAT THE VALUES OF  X(R)  ARE NON-DECREASING AND
C     DETERMINE
C     THE NUMBER  (MDIST)  OF DISTINCT VALUES OF  X(R).
C
      IERROR = 2
      MDIST = 1
      DO 40 R = 2, M
         IF (X(R).LT.X(R-1)) GO TO 380
         IF (X(R).EQ.X(R-1)) GO TO 40
         MDIST = MDIST + 1
   40 CONTINUE
C
C     IF THE  X(R)  ALL HAVE THE SAME VALUE, I.E.  MDIST = 1,
C     THE NORMALIZATION OF THE INDEPENDENT VARIABLE IS NOT
C     POSSIBLE.
C
      IERROR = 3
      IF (MDIST.EQ.1) GO TO 380
C
C     IF THE NUMBER OF DISTINCT VALUES OF  X(R)  FAILS TO EXCEED
C     THE MAXIMUM DEGREE  K  THERE IS NO UNIQUE POLYNOMIAL
C     APPROXIMATION OF THAT DEGREE.
C
      IERROR = 4
      IF (MDIST.LE.K) GO TO 380
C
C     CHECK THAT  NROWS  HAS BEEN SET SUFFICIENTLY LARGE.
C
      IERROR = 5
      IF (NROWS.LT.KPLUS1) GO TO 380
      IERROR = 0
C
      X1 = X(1)
      XM = X(M)
      D = XM - X1
C
C     THE INITIAL VALUES  WORK1(1, R)  (R = 1, 2, ..., M)  OF THE
C     WEIGHTED RESIDUALS AND THE VALUES  WORK1(2, R)  (R = 1, 2,
C     ..., M)
C     OF THE NORMALIZED INDEPENDENT VARIABLE ARE COMPUTED.  NOTE
C     THAT
C     WORK1(2, R)  IS COMPUTED FROM THE EXPRESSION BELOW RATHER
C     THAN THE
C     MORE NATURAL FORM
C
C     (2.0*X(R) - X1 - XM)/D,
C
C     SINCE THE FORMER GUARANTEES THE COMPUTED VALUE TO DIFFER FROM
C     THE TRUE VALUE BY AT MOST  4.0*MACHINE ACCURACY,  WHEREAS THE
C     LATTER HAS NO SUCH GUARANTEE.
C
      DO 60 R = 1, M
         WORK1(1,R) = W(R)*Y(R)
         WORK1(2,R) = ((X(R)-X1)-(XM-X(R)))/D
   60 CONTINUE
      I = 1
      BETAI = 0.0D0
      DO 360 IPLUS1 = 1, KPLUS1
C
C        SET STARTING VALUES FOR DEGREE  I.
C
         IPLUS2 = IPLUS1 + 1
         IF (IPLUS1.EQ.KPLUS1) GO TO 100
         DO 80 JPLUS1 = IPLUS2, KPLUS1
            A(IPLUS1,JPLUS1) = 0.0D0
   80    CONTINUE
         WORK2(1,IPLUS2) = 0.0D0
         WORK2(2,IPLUS2) = 0.0D0
  100    ALPIP1 = 0.0D0
         CI = 0.0D0
         DI = 0.0D0
         A(I,IPLUS1) = 0.0D0
         WORK2(1,IPLUS1) = 1.0D0
         IF (KPLUS1.GT.1) WORK2(2,1) = WORK2(1,2)
         DO 260 R = 1, M
            XCAPR = WORK1(2,R)
C
C           THE WEIGHTED VALUE  WORK1(3, R)  OF THE ORTHOGONAL
C           POLYNOMIAL OF
C           DEGREE  I  AT  X = X(R)  IS COMPUTED BY RECURRENCE FROM ITS
C           CHEBYSHEV-SERIES REPRESENTATION.
C
            IF (IPLUS1.GT.1) GO TO 120
            WRPR = W(R)*0.5D0*WORK2(1,1)
            WORK1(3,R) = WRPR
            GO TO 240
  120       J = IPLUS2
            IF (XCAPR.GT.0.5D0) GO TO 200
            IF (XCAPR.GE.-0.5D0) GO TO 160
C
C           GENTLEMAN*S MODIFIED RECURRENCE.
C
            FACTOR = 2.0D0*(1.0D0+XCAPR)
            DJ = 0.0D0
            BJ = 0.0D0
            DO 140 JREV = 1, I
               J = J - 1
               DJ = WORK2(1,J) - DJ + FACTOR*BJ
               BJ = DJ - BJ
  140       CONTINUE
            WRPR = W(R)*(0.5D0*WORK2(1,1)-DJ+0.5D0*FACTOR*BJ)
            WORK1(3,R) = WRPR
            GO TO 240
C
C           CLENSHAW*S ORIGINAL RECURRENCE.
C
  160       FACTOR = 2.0D0*XCAPR
            BJP1 = 0.0D0
            BJ = 0.0D0
            DO 180 JREV = 1, I
               J = J - 1
               BJP2 = BJP1
               BJP1 = BJ
               BJ = WORK2(1,J) - BJP2 + FACTOR*BJP1
  180       CONTINUE
            WRPR = W(R)*(0.5D0*WORK2(1,1)-BJP1+0.5D0*FACTOR*BJ)
            WORK1(3,R) = WRPR
            GO TO 240
C
C           REINSCH*S MODIFIED RECURRENCE.
C
  200       FACTOR = 2.0D0*(1.0D0-XCAPR)
            DJ = 0.0D0
            BJ = 0.0D0
            DO 220 JREV = 1, I
               J = J - 1
               DJ = WORK2(1,J) + DJ - FACTOR*BJ
               BJ = BJ + DJ
  220       CONTINUE
            WRPR = W(R)*(0.5D0*WORK2(1,1)+DJ-0.5D0*FACTOR*BJ)
            WORK1(3,R) = WRPR
C
C           THE COEFFICIENT  CI  OF THE  I TH  ORTHOGONAL POLYNOMIAL AND
C           THE
C           COEFFICIENTS  ALPIP1  AND  BETA I  IN THE
C           THREE-TERM RECURRENCE RELATION FOR THE ORTHOGONAL
C           POLYNOMIALS ARE COMPUTED.
C
  240       WRPRSQ = WRPR**2
            DI = DI + WRPRSQ
            CI = CI + WRPR*WORK1(1,R)
            ALPIP1 = ALPIP1 + WRPRSQ*XCAPR
  260    CONTINUE
         CI = CI/DI
         IF (IPLUS1.NE.1) BETAI = DI/DIM1
         ALPIP1 = 2.0D0*ALPIP1/DI
C
C        THE WEIGHTED RESIDUALS  WORK1(1, R)  (R = 1, 2, ..., M)  FOR
C        DEGREE  I  ARE COMPUTED, TOGETHER WITH THEIR SUM OF SQUARES,
C        SIGMAI.
C
         SIGMAI = 0.0D0
         DO 280 R = 1, M
            EPSR = WORK1(1,R) - CI*WORK1(3,R)
            WORK1(1,R) = EPSR
            SIGMAI = SIGMAI + EPSR**2
  280    CONTINUE
C
C        THE ROOT MEAN SQUARE RESIDUAL  S(I + 1)  FOR DEGREE  I  IS
C        THEORETICALLY UNDEFINED IF  M = I + 1  (THE CONDITION FOR THE
C        POLYNOMIAL TO PASS EXACTLY THROUGH THE DATA POINTS).  SHOULD
C        THIS
C        CASE ARISE THE R.M.S. RESIDUAL IS SET TO ZERO.
C
         IF (IPLUS1.GE.M) GO TO 300
         DF = M - IPLUS1
         S(IPLUS1) = SQRT(SIGMAI/DF)
         GO TO 320
  300    S(IPLUS1) = 0.0D0
C
C        THE CHEBYSHEV COEFFICIENTS  A(I + 1, 1), A(I + 1, 2), ...,
C        A(I + 1, I + 1)  TOGETHER WITH THE COEFFICIENTS
C        WORK2(1, 1), WORK2(1, 2), ..., WORK2(1, I + 1),   IN THE
C        CHEBYSHEV-SERIES REPRESENTATION OF THE  I TH  ORTHOGONAL
C        POLYNOMIAL ARE COMPUTED.
C
  320    DO 340 JPLUS1 = 1, IPLUS1
            JPLUS2 = JPLUS1 + 1
            PIJ = WORK2(1,JPLUS1)
            A(IPLUS1,JPLUS1) = A(I,JPLUS1) + CI*PIJ
            IF (JPLUS2.GT.KPLUS1) GO TO 380
            WORK2(1,JPLUS1) = WORK2(1,JPLUS2) + WORK2(2,JPLUS1) -
     *                        ALPIP1*PIJ - BETAI*WORK2(2,JPLUS2)
            WORK2(2,JPLUS2) = PIJ
  340    CONTINUE
         DIM1 = DI
         I = IPLUS1
  360 CONTINUE
  380 IF (IERROR) 400, 420, 400
  400 IFAIL = P01ABF(IFAIL,IERROR,SRNAME,0,P01REC)
      RETURN
  420 IFAIL = 0
      RETURN
      END
